#################################
# 12_estimateCJR_model_mexico.R #
#################################

rm(list=ls())

library (runjags)
library (coda)
library (reshape2)
library (random)
library (scales)

print (version)

# platform       aarch64-apple-darwin20      
# arch           aarch64                     
# os             darwin20                    
# system         aarch64, darwin20           
# status                                     
# major          4                           
# minor          1.2                         
# year           2021                        
# month          11                          
# day            01                          
# svn rev        81115                       
# language       R                           
# version.string R version 4.1.2 (2021-11-01)
# nickname       Bird Hippie     


# Set working directory
savePathRep <- "/all_data/"
setwd (savePathRep)

# load ("~/Dropbox/PartyEffects/JOP_replication/Step 4 results, figures and tables/fangyiRvectorJuly-corrected.RData")
# load("~/Dropbox/PartyEffects/MacPro/initialValues CJRvectorization20190710.Rda")
load ("mexico_rollcall_long_format_mean.Rdata")
load ("mexico_initial_values.Rda")

## Identification bills
potential.fix.bill <- c(22,89,144,150,174,217)


cjr.vector="model {
	for (i in 1:n.obs) {
		y[i] ~ dbern (pi[i])
		probit(pi[i]) <- beta[vote[i]]*theta[dep[i]]  + eta[vote[i]]*delta[dep[i]] - alpha[vote[i]]
	}
# PRIORS
for (j in 1:n.item){ alpha[j] ~ dnorm(0,1) }   
# Alternative identification: Beta (discrimination, dimension 1)
for(j in 1:(fix.bill[1]-1)) { beta[j]  ~ dnorm(0,1) }
beta[fix.bill[1]]~ dnorm(0,1)
for(j in (fix.bill[1]+1):(fix.bill[2]-1)) { beta[j]  ~ dnorm(0,1) }
beta[fix.bill[2]] <- 1.5
for(j in (fix.bill[2]+1):(fix.bill[3]-1)) { beta[j]  ~ dnorm(0,1) }
beta[fix.bill[3]] ~ dnorm(0,1)  #<- 2
for(j in (fix.bill[3]+1):(fix.bill[4]-1)) { beta[j]  ~ dnorm(0,1) }
beta[fix.bill[4]] <- -1
for(j in (fix.bill[4]+1):(fix.bill[5]-1)) { beta[j]  ~ dnorm(0,1) }
beta[fix.bill[5]] <-  2.5
for(j in (fix.bill[5]+1):(fix.bill[6]-1)) { beta[j]  ~ dnorm(0,1) }
beta[fix.bill[6]] ~ dnorm(0,1)
for(j in (fix.bill[6]+1):n.item) { beta[j]  ~ dnorm(0,1) }
# Eta (discrimination, dimension 2)
for(j in 1:(fix.bill[1]-1)) { eta[j]  ~ dnorm(0,1) }
eta[fix.bill[1]] ~ dnorm(0,1)
for(j in (fix.bill[1]+1):(fix.bill[2]-1)) { eta[j]  ~ dnorm(0,1) }
eta[fix.bill[2]] <- 2
for(j in (fix.bill[2]+1):(fix.bill[3]-1)) { eta[j]  ~ dnorm(0,1) }
eta[fix.bill[3]] ~ dnorm(0,1) #
for(j in (fix.bill[3]+1):(fix.bill[4]-1)) { eta[j]  ~ dnorm(0,1) }
eta[fix.bill[4]] <- -2 #
for(j in (fix.bill[4]+1):(fix.bill[5]-1)) { eta[j]  ~ dnorm(0,1) }
eta[fix.bill[5]] <- 0 #
for(j in (fix.bill[5]+1):(fix.bill[6]-1)) { eta[j]  ~ dnorm(0,1) }
eta[fix.bill[6]] ~ dnorm(0,1) #
for(j in (fix.bill[6]+1):n.item) { eta[j]  ~ dnorm(0,1) }
# ideal points
for(i in 1:n.legs)  { theta[i] ~ dnorm(0,1) }
for(i in 1:n.legs)  { delta[i] ~ dnorm(0,1) }
}"




cjr.data.vector <- dump.format(list(y=good.RC$RC
                                    , fix.bill=potential.fix.bill 
                                    , n.legs=max(good.comp.party$final.leg.no)
                                    , n.item=max(as.numeric(good.comp.party$vote))
                                    , n.obs=nrow(good.RC)
                                    , vote=as.numeric(good.comp.party$vote)
                                    , dep=good.comp.party$final.leg.no
))

cjr.parameters = c("theta","delta","alpha","beta","eta","deviance")


set.seed(2881)
eta1 = initials$eta + runif(length(initials$eta), -0.2, 0.2)
beta1 = initials$beta + runif(length(initials$beta), -0.2, 0.2)
eta1[potential.fix.bill[4:5]] <- c(NA,NA)
eta1[potential.fix.bill[2]] <- c(NA)
beta1[potential.fix.bill[2]] <- c(NA)
beta1[potential.fix.bill[4:5]] <- c(NA,NA)

set.seed(99292)
eta2 = initials$eta + runif(length(initials$eta), -0.2, 0.2)
beta2 = initials$beta + runif(length(initials$beta), -0.2, 0.2)

eta2[potential.fix.bill[4:5]] <- c(NA,NA)
eta2[potential.fix.bill[2]] <- c(NA)
beta2[potential.fix.bill[2]] <- c(NA)
beta2[potential.fix.bill[4:5]] <- c(NA,NA)


cjr.inits.1 <- function (){
  dump.format (
    list(
      theta = initials$theta + runif(length(initials$theta), -0.2, 0.2),
      delta = initials$delta + runif(length(initials$delta), -0.2, 0.2),
      alpha = initials$alpha + runif(length(initials$alpha), -0.2, 0.2),
      eta = eta1,
      beta = beta1
      ,'.RNG.name'="base::Wichmann-Hill"
      ,'.RNG.seed'=73828
    )
  )
}

cjr.inits.2 <- function (){
  dump.format (
    list(
      theta = initials$theta + runif(length(initials$theta), -0.2, 0.2),
      delta = initials$delta + runif(length(initials$delta), -0.2, 0.2),
      alpha = initials$alpha + runif(length(initials$alpha), -0.2, 0.2),
      eta = eta2,
      beta = beta2
      ,'.RNG.name'="base::Wichmann-Hill"
      ,'.RNG.seed'=88393
    )
  )
}

cjr.model <- run.jags(
  model=cjr.vector,
  monitor=cjr.parameters,
  method="parallel",
  n.chains=2,
  data=cjr.data.vector,
  inits=list (cjr.inits.1(), cjr.inits.2()),
  thin=50, burnin=5000, sample=100,
  check.conv=FALSE, plots=FALSE)

chainsCJR <- mcmc.list(list (cjr.model$mcmc[[1]], cjr.model$mcmc[[2]]))  
gelman.CJR <- gelman.diag (chainsCJR, multivariate=F) 

# Uncomment next line to save results
# save(chainsCJR, file="mexico_CJR_auxiliary_output.Rda")
